function [sgb, ...
          equalorder, ...
          Ugridreduced, sgreduced,...
          id_Ugridreduced,...
          s_id_gridreduced, sU_id_gridreduced,...
          n_U_sets, Tcoord, indices_pairs,...
          numeqconstraints,numzeros,coordrowseq, fillAeq]=useful_anyparam2(u0grid, u1grid, u2grid, sg)     
%% U sets
n_U_sets=3; %number of U-sets (dimension of the distribution of the taste shock differences)
Ugrid=cell(1,n_U_sets); %cell 1xn_U_sets

Ugrid{1}=[u0grid-u1grid -u0grid+u1grid zeros(sg,1) Inf(sg,1) -Inf(sg,1)];
Ugrid{2}=[u0grid-u2grid -u0grid+u2grid zeros(sg,1) Inf(sg,1) -Inf(sg,1)]; 
Ugrid{3}=[Inf(sg,1) -Inf(sg,1) u2grid-u1grid -u2grid+u1grid zeros(sg,1)];

sgb=0;     
%% Reduce number of grid points to check by applying our partitioning arguments
%equalorder (sgx1) assigns the same index to parameters vectors such that Ugrid{1}, Ugrid{2}, ..., Ugrid{n_U_sets} have the same order

%%%%%%%% Criterion 1 %%%%%%%%
d=cell(n_U_sets,1);
index=cell(n_U_sets,1);
for j=1:n_U_sets 
    [sortedUjgrid,index{j}] = sort(Ugrid{j},2); 
    %sortedUjgrid gives Ugrid{j} with each row sorted from smallest to largest; it is sgxn_columns_Ujgrid
    %for each row i of sortedUjgrid and for k=1,2,...,n_columns_Ujgrid, index{j}(i,k) gives the column position in Ugrid{j}(i,:) of sortedUjgrid(i,k); it is sgxn_columns_Ujgrid
    d{j} = diff(sortedUjgrid,[],2) == 0; %sg x (n_columns_Ujgrid-1)
    %Suppose that Ugrid{j} has 3 columns. Then, for each row i of sortedUjgrid, d{j} will be: 
    %[1 1] if sortedUgridj(i,2)-sortedUgridj(i,1)=0 and sortedUgridj(i,3)-sortedUgridj(i,2)=0
    %[1 0] if sortedUgridj(i,2)-sortedUgridj(i,1)=0 and sortedUgridj(i,3)-sortedUgridj(i,2)~=0
    %[0 1] if sortedUgridj(i,2)-sortedUgridj(i,1)~=0 and sortedUgridj(i,3)-sortedUgridj(i,2)=0
    %[0 0] if sortedUgridj(i,2)-sortedUgridj(i,1)~=0 and sortedUgridj(i,3)-sortedUgridj(i,2)~=0
end
%%%%%%%% Criterion 2 %%%%%%%%
n_columns_Ujgrid=5;
Ugrid_sum=zeros(sg,n_columns_Ujgrid+n_columns_Ujgrid^2);
for g=1:sg
Ugrid_sum(g,:)=[Ugrid{1}(g,:) reshape((Ugrid{2}(g,:)+Ugrid{3}(g,:).')',[],1)'];
end

[sortedUgrid_sum,index_sum] = sort(Ugrid_sum,2); 
d_sum = diff(sortedUgrid_sum,[],2) == 0; %sg x (n_columns_Ujgrid*n_U_sets-1)

%%%%%%%% Combine %%%%%%%%
[~,~,equalorder] = unique([index_sum d_sum [index{:}] [d{:}] ],'rows', 'stable'); %sgx1 

%Reduced grid of parameter vectors and U sets
[~,index_final,~]=unique(equalorder, 'stable'); %one row from equalorder per ordering group together with its index
Ugridreduced=cell(1,n_U_sets);
for j=1:n_U_sets 
    Ugridreduced{j}=Ugrid{j}(index_final,:); %sgreduced x n_columns_Ujgrid
end
sgreduced=size(Ugridreduced{j},1); %size of the reduced grid of parameter values


%% Construct indices to keep track of equal elements inside Ugridreduced{j} for j=1,...,n_U_sets
id_Ugridreduced=cell(n_U_sets,1);
for j=1:n_U_sets 
    id_Ugridreduced{j} = NaN(size(Ugridreduced{j})); % %sgreduced x n_columns_Ujgrid
    for k = 1:sgreduced
        [~, ~, id_Ugridreduced{j}(k,:)] = unique(Ugridreduced{j}(k,:), 'stable'); 
    end
end

s_id_gridreduced=zeros(sgreduced, n_U_sets); %useful to set number of unknowns of the linear programming
for j=1:n_U_sets
    s_id_gridreduced(:,j)=max(id_Ugridreduced{j},[],2); %number of distinct elements of Ugridreduced{j}
end
sU_id_gridreduced=prod(s_id_gridreduced,2); %sgreducedx1 reporting the cardinality of the Cartesian product of the U sets
                                            %sU_id_gridreduced(g) is the number of unknowns of the LP for each g-th parameter value
s_id_gridreduced_2=ones(1, n_U_sets)*n_columns_Ujgrid; %number of elements of Ugridreduced{j} (count as separate even equal elements)
sU_id_gridreduced_2=prod(s_id_gridreduced_2); %cardinality of the Cartesian product of the U sets     
  

%% Invariant components of the linear programming problem 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% Inequality constraints %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

vectors = arrayfun(@(x) {1:x}, s_id_gridreduced_2); 
n = numel(vectors);
Tcoord_temp = cell(1,n);
[Tcoord_temp{:}] = ndgrid(vectors{:}); 
Tcoord_temp = cat(n+1, Tcoord_temp{:});
Tcoord = reshape(Tcoord_temp,[],n); %matrix of dimension sU_id_gridreduced_2xn_U_sets reporting the column coordinates of all the possible n_U_sets-tuples from the U-sets
%For example if n_U_sets=3 and s_id_gridreduced_2 is [2 3 1], then 
%[ca, cb, cc] = ndgrid(1:2, 1:3, 1:1);
%Tcoord = [ca(:), cb(:), cc(:)];   
%Tcoord= [1     1     1
%         2     1     1
%         1     2     1
%         2     2     1
%         1     3     1
%         2     3     1];
indices_pairs=pairIndices(sU_id_gridreduced_2); %row indices of all possible unordered pairs of rows from the matrix obtained by taking the Cartesian product of the U sets
                                         %[sU_id_gridreduced_2*(sU_id_gridreduced_2-1)/2]x[2]   
%Continuing the example above
%indices_pairs=[1 1; 1 2; ...; 1 27; 2 3; ...; 2 27; ... 26 27]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% Equality constraints %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

numeqconstraints=3+(5+2*sgb+5+2*sgb-1)*(5+2*sgb-1)+(5+2*sgb)^2+1+(2+sgb)*3;  
numzeros=(5+2*sgb+5+2*sgb-1)*(5+2*sgb-1)+(5+2*sgb)^2;


coordrowseq=[1 1 1 ... %observational equivalence
             2 2   ...
             3     ...
             4:1:3+numzeros ... %zeros
             3+numzeros+1 ... %ones
             kron((3+numzeros+1+1:1:3+numzeros+1+(1+sgb)*3), ones(1,2)) ... %symmetry
             3+numzeros+1+(1+sgb)*3+1:1:3+numzeros+1+(1+sgb)*3+3 ];  

                         
fillAeq=[1 -1 -1 ... %observational equivalence
         1 -1    ...
         1       ...
         ones(1,numzeros) ... %zeros
         1 ...%ones 
         repmat([1 1],1, (1+sgb)*3) ...%symmetry
         ones(1,3)];    
end